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Abstract: Extending the Liouville-Caputo definition of a fractional derivative to a nonlocal 
covariant generalization of arbitrary bound operators acting on multidimensional Riemannian 
spaces an appropriate approach for the 3D shape recovery of aperture afflicted 2D slide sequences 
is proposed. We demonstrate, that the step from a local to a nonlocal algorithm yields an order 
of magnitude in accuracy and by using the specific fractional approach an additional factor 2 
in accuracy of the derived results. 

Keywords: Fractional calculus, computer graphics, image processing, shape recovery, confocal 
microscopy, modified Laplacian. 



1. INTRODUCTION 

Prom a historical point of view, fractional calculus provides 

us with a set of axioms and methods to extend the concept 
of a derivative operator from integer order n to arbitrary 
order a, where a is a real or complex value. 
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(1) 



In the sense of (1) fractional calculus has been frequently 
applied in the area of image processing, see e.g. [Falzon 
(1994)], [Ortigueira (2003)], [Sparavigna (2009)]. 

Alternatively we may consider fractional calculus as a 
specific prescription to extend the definition of a local 
operator to the nonlocal case. In this lecture, we will 
present a covariant, multidimensional generalization of the 
fractional derivative definition, which may be applied to 
any bound operator on the Riemannian space. 

As a first application, we will propose a specific non-local 
extension of the modified local Laplace-operator, which 
is widely used in problems of image processing. We will 
especially compare the local to the nonlocal approach for 
3D-shape recovery from a set of 2D aperture afflicted 
slide sequences, which may be obtained e.g. in confo- 
cal microscopy or autofocus algorithms [Zernike (1935)], 
[Spencer (1982)]. 

It will be shown, that a major improvement of results is 
achieved for the nonlocal version of the modified Laplace- 
operator. 

2. THE GENERALIZED FRACTIONAL DERIVATIVE 



of arbitrary local operators. For that purpose, we start 
with the Liouville definition of the left and right fractional 
integral [Liouville (1832)]: 
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With a slight modification of the fractional parameter 
a = 1 — a, where a is in the interval < a < 1. 
Consequently for the limiting case a = 1+ and /_ both 
coincide with the \mit-operator and for a = I and 7_ 
both correspond to the standard integral operator. 

/+ and /_ may be combined to define a regularized 
Liouville integral [Herrmann (2011)]: 
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where we have introduced the symmetric shift-operator: 
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We will propose a reinterpretation of the fractional cal- 
culus as a specific procedure for a non-local extension The regularized fractional Liouville-Caputo derivative may 
7 now be defined as: 



d"^ f{x) = I°'d.J{x) (7) /"(a;,?/) explicitly reads: 



with the abbreviation dxf{x) = / {x). This definition of a 
fractional derivative coincides with Feller's [Feller (1952)] 
definition pd^iO) for the special case 6 — 1. 

We may interpret /" as a non-localization operator, which 
is applied to the local derivative operator to determine a 
specific non-local extension of the same operator. There- 
fore the fractional extension of the derivative operator is 
separated into a sequential application of the standard 
derivative followed by a non-localization operation. The 
classical interpretation of a fractional integral is changed 
from the inverse operation of a fractional derivative to a 
more general interpretation of a non-localization proce- 
dure, which may be easily interpreted in the area of image 
processing as a blur effect. 

This is a conceptual new approach, since it may be easily 
extended to other operators e.g. higher order derivatives 
or space dependent operators, e.g. for we obtain: 

{dirf{x)=i"dif{x) (11) 

^^l^^d^e-'m)dlf{x) (12) 
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(15) 

which is nothing else but the Riesz [Riesz (1949)] definition 

of a fractional derivative. 

Therefore we define the following fractional extension of a 
local operator locaiO to the non-local case 

noniocaiO" fix) = /"localO/(a;) (16) 

as the covariant generalization of the Liouvillc-Caputo 
fractional derivative to arbitrary operators on TZ. 

This definition may be easily extended to the multidimen- 
sional case, interpreting the variable ^ as a measure of 

distance. 



In two dimensions, with 



and with 
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< a < 2 



(20) 



which is normalized such, that the eigenvalue spectrum 
for: 

I"f{x,y) = Kfix,y) (21) 
with the eigenfunctions f{x, y) = exp''^^^+'*'^^ follows as: 



K=(fc2 + fc2)-a/2 



(22) 



It should be noted, that the validity range for a spans from 

< a < 2, since we deal with a two-dimensional problem. 
Obviously within the framework of signal processing, the 
non-localization operator may be interpreted as a low-pass 
filter. 

In the following sections, we will use this operator for a well 

defined extension of the standard algorithm used for 3D- 
shape recovery from aperture afflicted 2D-slide sequences 
to a generalized, fractional nonlocal version, which results 
in a very stable procedure with drastically reduced errors. 

We will first present the minimal standard method and its 
limitations in the next section. 

3. THE LOCAL APPROACH 

In a set {pi{zi), i = 1, N} of N 2D-slides with increasing 
focal distance Zi, {zi,\/zi : z-^^^ < Zi < Zmin,* = l,...,iV} 
every slide contains areas with focused as well as defocused 
parts of the specimen considered. In the first row of Fig. 

1 we present two examples from a slide-sequence of a 
spherical object with radius r = 1 located at z = in the 
x,y-plane, where the focal plane was chosen to be z = 0.45 
and z = 0.8 respectively. 

For a 3D-shape recovery in a first step for a given slide the 
parts being in focus have to be extracted. For a textured 
object, areas in focus are dominated by a larger amount of 

high frequency contributions, while for out of focus parts 
mainly the low frequency amount of texture survives. 

Consequently an appropriate operator to determine the 
high-frequency domains is the modified Laplacian locai^ 
given e.g. by: 

(i_iA/)(a;,y) = (|J| + |J|)/(a.,y) (23) 

X y 

where 1 1 denotes the absolute value. 

In the discrete case with a symmetrically discretized 
function fij on a rectangular domain aJmin < a; < Xmax 
and t/min <y<y max- 

fij = f{Xmin + -ih, ymin + jh) (24) 
i ~ 0) •••) *max ) J ~ 0) •••) Jmax 

with stepsize h in both x- and y-direction, the same 
operator is given by: 




nonlocal ^ *^ 

z=0.45 z=0.80 

Fig. 1. Application of the local and nonlocal modified 
Laplace-operator to 2 different original slides (left 
column z — 0.45 and right column z — 0.8) from 
a 2D-slide sequence z g {0,1}. From top to bottom 
original slide, result of local modified Laplacian from 
(25) iocaiA((7 = 1) and result of nonlocal operator are 
shown. 

(localA(g) f)ij = 

I fi+q,j ~ ^/'' j' + fi-q,j I I I fi,j+q j" + fi-J-g i 

* = •••7«.„ax -Q, j = q, ■•■,imax - g (25) 

and elsewhere, where the free parameter q has to 
be chosen according to the Nyquist-Shannon sampling 
theorem [Shannon (1949)] to be of order of the inverse 
average wavelength cj of the texture applied to the object 
considered 

qh^2/uj (26) 

a requirement, which can be fulfilled only locally for 
random generated textures and for regular textures on 
curved surfaces respectively. 

An application of the modified Laplacian to every slide in 
a set {pkizk)} leads to a set of intensity values {pij{q, z^.)} 
at a given pixel-position at ij: 

loc^,lM<l){Pi3.ki^k)} = {pijiq,Zk)} (27) 
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Fig. 2. Comparison of recovered positions with original 
height positions for a set of points along the y-axis. 
Black circles mark the correct positions, lines repre- 
sent recovered positions based on the local modified 
Laplacian iocaiA((7) from (25) for different values of 
step-size q. It should be noted, that there is no fixed 
value of q, which uniquely may be used to determine 
all positions. There are drop outs for every curve. 
Errors are listed in Table 1. 

In the second row of Fig. 1 the result of an application of 
the discrete modified Laplacian with g = 1 to the original 
slides presented in the first row, is demonstrated. 

It is assumed, that for a fixed q a maximum exists in 
Pij{Q,Zk) for a given k. A parabolic fit of {pij{q, z^)} near 
k helps to determine the position optz{q), where pij{q,z) 
is maximal: 

""^ 2 p^j{q, Zf^^^) - 2pij{q, z^) + pij{q, z^_^) 

(28) 

In the center row of Fig. 1 we present the result of the 
application of (25) onto the original slides. The gray-level 
indicates the intensity values {pij(q, Zk)} for Zk — 0.45 and 
Zfc = 0.8 respectively. In Fig 2 recovered Zopt('?) along the 
positive y-axis are compared for different q- values with the 
original height-values. 

Obviously there is no unique optimum choice for q, which 
works for all positions simultaneously. The proposed sim- 
ple local approach is not very effective, instead it generates 
drop outs as a result of an interference of varying texture 
scaling with the fixed step size q. For a realistic treatment 
of 3D-shape recovery a more sophisticated procedure is 
necessary. 

Consequently a nonlocal approach, which weights the 
different contributions for a varying step size is a promising 
and well defined approach. Indeed it will enhance the 
quality of the results significantly, as will be demonstrated 
in the next section. 

4. THE NONLOCAL APPROACH 

The generalized fractional approach extends the above pre- 
sented local algorithm. The nonlocal modified Laplacian 
according to (16) is given by : 



Table 1. Comparison of rms-errors in % for 
nonlocal modified Laplacian from (30) for dif- 
ferent a to local approach from (27) with vary- 
ing q in the last column. 
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0.277 


0.265 
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0.215 
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4.228 
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(29) 

Therefore we obtain a well defined two step procedure. 
First, the local operator is applied, followed by the nonlo- 
calization integral. 

In the discrete case, applying the nonlocal Laplacian to 
every slide in a slide set, the first step is therefore identical 
with (27) and yields a set of intensity values {pij{q, z^)} at 
a given pixel-position at ij. An application of the discrete 
version of then leads to: 

nonlocalA"(9)pij,fc(2fe) = /"local A ((j) Kj^fe (Zfc ) 
1 

~ 2«-2r(a/2)2sin(a7r/2) ^ (30) 

/ da / d6(C?+^2)^^""'^5(a,6)p.,(g,^fe) 
Jo Jo 

= p?jiq,zk) 

where we have introduced a cutoff (, which limits the 
integral on the finite domain of pixel values. If we interpret 
the intensity values as constant function values at position 
ij with size h, the integration may be performed fully 
analytically. In the appendix we have listed the resulting 
matrix-operator for C = 4. 

The resulting nonlocal intensities p^ (g, z^,) are presented 
in the lower row of Fig. 1. The nonlocal approach reduces 
the granularity of the local operator and a more smooth 
behaviour of intensities results. 

Since this is the only modification of the local approach, 
the recovery of the height information for every pixel is 
similar to (28) 

op.na, q) 2 pg(g, z-,^^) - 2pf^.(g, z^) + pf^.(g, J 

(31) 

In Fig. 3 results are plotted for different values of a. 
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Fig. 3. Comparison of recovered positions with original 
height positions for a set of points along the y- 
axis. Black circles mark the correct positions, lines 
represent recovered positions based on the nonlocal 
modified Laplacian iocaiA((7) from (30) for different 
values of the fractional parameter a. The algorithm is 
very stable against a variation of a. In the limit a = 
the nonlocal approach reduces to the local scenario. 
Errors are listed in Table 1. 



In Table 1 a listing of errors is given for the local and the 
nonlocal algorithm presented. 

We may conclude, that the nonlocal approach is very 
robust and stable in a wide range of a and C, values 
respectively. We gain one order of magnitude in accuracy 
using the nonlocal modified Laplacian. An additional 
factor 2 in accuracy is obtained if we chose the optimal 
fractional {a, C} parameter set. 



5. APPENDIX 

The discrete version of the nonlocalization operator 7" 
from (30) may be interpreted as a matrix operation M{a) 
on fif 

= Af(a)/„- (32) 

M{a) is a quadratic (2^ 1) x (2^ + 1) matrix with the 
symmetry properties 

M(a),j = Af(a)_,,j = M(a)_,,_j = M{a),^_j = ^//(a)^-, 

(33) 

Setting the normalization condition M{q)oq = 1 the 
integral may be solved analytically for stepwise constant 
pixel values pij. As an example, we present the fourth 
quadrant of M{a) for ^ = 4 in units h: 



M(2) : 



M(3/2) : 



M(l): 



/I. 1. 1. 1. l.N^ 

1. 1. 1. 1. 1 

1. 1. 1. 1. 1 

1. 1. 1. 1. 1, 

Vi. 1. 1. 1. 1./ 

/ 1. 0.570351 
0.570351 0.478687 
0.400990 0.379117 
0.326971 0.318443 

Vo. 283027 0.278761 

/ 1. 0.294441 
0.294441 0.205559 
0.143268 0.127951 
0.094982 0.090073 

\0.071095 0.068963 



(34) 



0.400990 0.326971 0.283027\ 
0.379117 0.318443 0.278761 
0.336822 0.298160 0.267640 
0.298160 0.274801 0.253092 
0.267640 0.253092 0.237922/ 

0.143268 0.094982 0.071095\ 
0.127951 0.090073 0.068963 
0.100830 0.078927 0.063559 
0.078927 0.067014 0.056825 
0.063559 0.056825 0.050208/ 



Spencer, M. (1982). Fundamentals of Light Microscopy, 

Cambridge University Press. 
Zernike, M. (1935). Das Phasenkontrastverfaren bei der 

mikroskopischen Beobachtung Z. Tech. Phys. 16, 454- 

457. 



M(l/2) : 



M(0) : 



^ 1. 0.116147 0.038486 0.020685 0.013374\ 

0.116147 0.066866 0.032440 0.019096 0.012776 ' 

0.038486 0.032440 0.022637 0.015652 0.011301 

0.020685 0.019096 0.015652 0.012237 0.009550 

Vo.013374 0.012776 0.011301 0.009550 0.007929/ 

/I. 0. 0. 



0. 
0. 
0. 

Vo. 



0. 

0. 

0. 

0. O.J 



(35) 



Obviously there is a smooth transition from a local 
(q: = 0) to a more and more nonlocal opc;ration, which in 
the limiting case (a = 2) may be interpreted as the result 
of the use of a pinhole camera with finite hole radius 
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